# Replication code for Claudia J. Kim and Taylor C. Boas, "Activist Disconnect: Social Movements, Public Opinion, and U.S. Military Bases in East Asia," Armed Forces and Society.

# Analysis conducted in R 3.6.0 on MacOS 10.13.6

# NOTE: This file reproduces Appendix Figures 9-10. 

# Set working directory as appropriate
# setwd('~/Dropbox/us_bases_experiment/replication/')

# Clean desktop and load packages. Please make sure all necessary packages are installed.

rm(list=ls(all=T))

library(Hmisc)

# Load cleaned dataset

load('basedata_cleaned.RData')

# Function to demean covariates

demean<-function(x) x-mean(x,na.rm=T)

# Rename/recode variables

basedata$Age<-basedata$Q4
basedata$Male<-basedata$Q33==1
basedata$Education<-basedata$Q34
basedata$Satisfaction<-basedata$Q35
basedata$Duty_vote<-basedata$Q6
basedata$Efficacy<-basedata$Q7
basedata$Petition<-basedata$Q8
basedata$Boycott<-basedata$Q9
basedata$Demonstration<-basedata$Q10
basedata$Donate<-basedata$Q11
basedata$Ideology<-basedata$Q13
basedata$Extremism<-abs(basedata$Q13-5.5)
basedata$Distance_org<-abs(basedata$Q13-ifelse(basedata$Country=='Korea',basedata$Q16,basedata$Q20))
basedata$Local_news<-basedata$Q22
basedata$Nat_news<-basedata$Q23
basedata$Conflict<-basedata$Q28
basedata$Intervention<-basedata$Q29
basedata$Attitude<-basedata$Q30
basedata$Behavior<-basedata$Q37
basedata$TreatA<-basedata$Q30_cond==1
basedata$TreatB<-basedata$Q30_cond==2
basedata$TreatC<-basedata$Q30_cond==3
basedata$TreatD<-basedata$Q30_cond==4
basedata$TreatE<-basedata$Q30_cond==5
basedata$TreatF<-basedata$Q30_cond==6
basedata$TreatG<-basedata$Q30_cond==7
basedata$TreatH<-basedata$Q30_cond==8
basedata$Daegu<-basedata$Q5=='Daegu'
basedata$Gyeonggi<-basedata$Q5=='Gyeonggi'
basedata$Kanagawa<-basedata$Q5=='Kanagawa'
basedata$Okinawa<-basedata$Q5=='Okinawa'
basedata$Region<-basedata$Q5
basedata$Screener<-basedata$Q25==14
basedata$Mil_presence<- basedata$Q32

# Fixed effects and covariates used as controls in regressions

regions<-c('demean(Daegu)', 'demean(Gyeonggi)', 'demean(Kanagawa)', 'demean(Okinawa)')

covars<-c('demean(Age)', 'demean(Male)', 'demean(Education)', 'demean(Satisfaction)', 'demean(Duty_vote)', 'demean(Efficacy)', 'demean(Petition)', 'demean(Boycott)', 'demean(Demonstration)', 'demean(Donate)', 'demean(Ideology)', 'demean(Extremism)', 'demean(Distance_org)', 'demean(Local_news)', 'demean(Nat_news)', 'demean(Conflict)', 'demean(Intervention)')

############
# Testing H1
############

# No Covariates

h1_pure_att_screen<-with(basedata[basedata$Screener & (basedata$TreatC | basedata$TreatA),], lm(as.formula(paste0('Attitude ~ TreatC + ', paste(regions, '* TreatC', collapse = ' + ')))))
h1_reg_att_screen<-with(basedata[basedata$Screener & (basedata$TreatC | basedata$TreatB),], lm(as.formula(paste0('Attitude ~ TreatC + ', paste(regions, '* TreatC', collapse = ' + ')))))
h1_pure_beh_screen<-with(basedata[basedata$Screener & (basedata$TreatC | basedata$TreatA),], lm(as.formula(paste0('Behavior ~ TreatC + ', paste(regions, '* TreatC', collapse = ' + ')))))
h1_reg_beh_screen<-with(basedata[basedata$Screener & (basedata$TreatC | basedata$TreatB),], lm(as.formula(paste0('Behavior ~ TreatC + ', paste(regions, '* TreatC', collapse = ' + ')))))

# Covariates

h1_pure_att_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatC | basedata$TreatA),], lm(as.formula(paste0('Attitude ~ TreatC + ', paste(covars, '* TreatC', collapse = ' + '), '+', paste(regions, '* TreatC', collapse = ' + ')))))
h1_reg_att_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatC | basedata$TreatB),], lm(as.formula(paste0('Attitude ~ TreatC + ', paste(covars, '* TreatC', collapse = ' + '), '+', paste(regions, '* TreatC', collapse = ' + ')))))
h1_pure_beh_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatC | basedata$TreatA),], lm(as.formula(paste0('Behavior ~ TreatC + ', paste(covars, '* TreatC', collapse = ' + '), '+', paste(regions, '* TreatC', collapse = ' + ')))))
h1_reg_beh_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatC | basedata$TreatB),], lm(as.formula(paste0('Behavior ~ TreatC + ', paste(covars, '* TreatC', collapse = ' + '), '+', paste(regions, '* TreatC', collapse = ' + ')))))

############
# Testing H2
############

# No Covariates

h2_pure_att_neg_screen<-with(basedata[basedata$Screener & (basedata$TreatD | basedata$TreatA) & basedata$Region!='Daegu',], lm(as.formula(paste0('Attitude ~ TreatD + ', paste(regions, '* TreatD', collapse = ' + ')))))
h2_reg_att_neg_screen<-with(basedata[basedata$Screener & (basedata$TreatD | basedata$TreatB) & basedata$Region!='Daegu',], lm(as.formula(paste0('Attitude ~ TreatD + ', paste(regions, '* TreatD', collapse = ' + ')))))
h2_pure_beh_neg_screen<-with(basedata[basedata$Screener & (basedata$TreatD | basedata$TreatA) & basedata$Region!='Daegu',], lm(as.formula(paste0('Behavior ~ TreatD + ', paste(regions, '* TreatD', collapse = ' + ')))))
h2_reg_beh_neg_screen<-with(basedata[basedata$Screener & (basedata$TreatD | basedata$TreatB) & basedata$Region!='Daegu',], lm(as.formula(paste0('Behavior ~ TreatD + ', paste(regions, '* TreatD', collapse = ' + ')))))

h2_pure_att_pos_screen<-with(basedata[basedata$Screener & (basedata$TreatD | basedata$TreatA) & basedata$Region=='Daegu',], lm(as.formula(paste0('Attitude ~ TreatD + ', paste(regions, '* TreatD', collapse = ' + ')))))
h2_reg_att_pos_screen<-with(basedata[basedata$Screener & (basedata$TreatD | basedata$TreatB) & basedata$Region=='Daegu',], lm(as.formula(paste0('Attitude ~ TreatD + ', paste(regions, '* TreatD', collapse = ' + ')))))
h2_pure_beh_pos_screen<-with(basedata[basedata$Screener & (basedata$TreatD | basedata$TreatA) & basedata$Region=='Daegu',], lm(as.formula(paste0('Behavior ~ TreatD + ', paste(regions, '* TreatD', collapse = ' + ')))))
h2_reg_beh_pos_screen<-with(basedata[basedata$Screener & (basedata$TreatD | basedata$TreatB) & basedata$Region=='Daegu',], lm(as.formula(paste0('Behavior ~ TreatD + ', paste(regions, '* TreatD', collapse = ' + ')))))

# Covariates

h2_pure_att_neg_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatD | basedata$TreatA) & basedata$Region!='Daegu',], lm(as.formula(paste0('Attitude ~ TreatD + ', paste(covars, '* TreatD', collapse = ' + '), '+', paste(regions, '* TreatD', collapse = ' + ')))))
h2_reg_att_neg_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatD | basedata$TreatB) & basedata$Region!='Daegu',], lm(as.formula(paste0('Attitude ~ TreatD + ', paste(covars, '* TreatD', collapse = ' + '), '+', paste(regions, '* TreatD', collapse = ' + ')))))
h2_pure_beh_neg_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatD | basedata$TreatA) & basedata$Region!='Daegu',], lm(as.formula(paste0('Behavior ~ TreatD + ', paste(covars, '* TreatD', collapse = ' + '), '+', paste(regions, '* TreatD', collapse = ' + ')))))
h2_reg_beh_neg_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatD | basedata$TreatB) & basedata$Region!='Daegu',], lm(as.formula(paste0('Behavior ~ TreatD + ', paste(covars, '* TreatD', collapse = ' + '), '+', paste(regions, '* TreatD', collapse = ' + ')))))

h2_pure_att_pos_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatD | basedata$TreatA) & basedata$Region=='Daegu',], lm(as.formula(paste0('Attitude ~ TreatD + ', paste(covars, '* TreatD', collapse = ' + '), '+', paste(regions, '* TreatD', collapse = ' + ')))))
h2_reg_att_pos_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatD | basedata$TreatB) & basedata$Region=='Daegu',], lm(as.formula(paste0('Attitude ~ TreatD + ', paste(covars, '* TreatD', collapse = ' + '), '+', paste(regions, '* TreatD', collapse = ' + ')))))
h2_pure_beh_pos_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatD | basedata$TreatA) & basedata$Region=='Daegu',], lm(as.formula(paste0('Behavior ~ TreatD + ', paste(covars, '* TreatD', collapse = ' + '), '+', paste(regions, '* TreatD', collapse = ' + ')))))
h2_reg_beh_pos_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatD | basedata$TreatB) & basedata$Region=='Daegu',], lm(as.formula(paste0('Behavior ~ TreatD + ', paste(covars, '* TreatD', collapse = ' + '), '+', paste(regions, '* TreatD', collapse = ' + ')))))

########################
# Testing H3a--pragmatic
########################

# No Covariates

h3a_prag_pure_att_screen<-with(basedata[basedata$Screener & (basedata$TreatE | basedata$TreatA),], lm(as.formula(paste0('Attitude ~ TreatE + ', paste(regions, '* TreatE', collapse = ' + ')))))
h3a_prag_reg_att_screen<-with(basedata[basedata$Screener & (basedata$TreatE | basedata$TreatB),], lm(as.formula(paste0('Attitude ~ TreatE + ', paste(regions, '* TreatE', collapse = ' + ')))))
h3a_prag_pure_beh_screen<-with(basedata[basedata$Screener & (basedata$TreatE | basedata$TreatA),], lm(as.formula(paste0('Behavior ~ TreatE + ', paste(regions, '* TreatE', collapse = ' + ')))))
h3a_prag_reg_beh_screen<-with(basedata[basedata$Screener & (basedata$TreatE | basedata$TreatB),], lm(as.formula(paste0('Behavior ~ TreatE + ', paste(regions, '* TreatE', collapse = ' + ')))))

# Covariates

h3a_prag_pure_att_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatE | basedata$TreatA),], lm(as.formula(paste0('Attitude ~ TreatE + ', paste(covars, '* TreatE', collapse = ' + '), '+', paste(regions, '* TreatE', collapse = ' + ')))))
h3a_prag_reg_att_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatE | basedata$TreatB),], lm(as.formula(paste0('Attitude ~ TreatE + ', paste(covars, '* TreatE', collapse = ' + '), '+', paste(regions, '* TreatE', collapse = ' + ')))))
h3a_prag_pure_beh_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatE | basedata$TreatA),], lm(as.formula(paste0('Behavior ~ TreatE + ', paste(covars, '* TreatE', collapse = ' + '), '+', paste(regions, '* TreatE', collapse = ' + ')))))
h3a_prag_reg_beh_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatE | basedata$TreatB),], lm(as.formula(paste0('Behavior ~ TreatE + ', paste(covars, '* TreatE', collapse = ' + '), '+', paste(regions, '* TreatE', collapse = ' + ')))))


##########################
# Testing H3a--ideological
##########################

# No Covariates

h3a_ideol_pure_att_screen<-with(basedata[basedata$Screener & (basedata$TreatF | basedata$TreatA),], lm(as.formula(paste0('Attitude ~ TreatF + ', paste(regions, '* TreatF', collapse = ' + ')))))
h3a_ideol_reg_att_screen<-with(basedata[basedata$Screener & (basedata$TreatF | basedata$TreatB),], lm(as.formula(paste0('Attitude ~ TreatF + ', paste(regions, '* TreatF', collapse = ' + ')))))
h3a_ideol_pure_beh_screen<-with(basedata[basedata$Screener & (basedata$TreatF | basedata$TreatA),], lm(as.formula(paste0('Behavior ~ TreatF + ', paste(regions, '* TreatF', collapse = ' + ')))))
h3a_ideol_reg_beh_screen<-with(basedata[basedata$Screener & (basedata$TreatF | basedata$TreatB),], lm(as.formula(paste0('Behavior ~ TreatF + ', paste(regions, '* TreatF', collapse = ' + ')))))

# Covariates

h3a_ideol_pure_att_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatF | basedata$TreatA),], lm(as.formula(paste0('Attitude ~ TreatF + ', paste(covars, '* TreatF', collapse = ' + '), '+', paste(regions, '* TreatF', collapse = ' + ')))))
h3a_ideol_reg_att_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatF | basedata$TreatB),], lm(as.formula(paste0('Attitude ~ TreatF + ', paste(covars, '* TreatF', collapse = ' + '), '+', paste(regions, '* TreatF', collapse = ' + ')))))
h3a_ideol_pure_beh_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatF | basedata$TreatA),], lm(as.formula(paste0('Behavior ~ TreatF + ', paste(covars, '* TreatF', collapse = ' + '), '+', paste(regions, '* TreatF', collapse = ' + ')))))
h3a_ideol_reg_beh_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatF | basedata$TreatB),], lm(as.formula(paste0('Behavior ~ TreatF + ', paste(covars, '* TreatF', collapse = ' + '), '+', paste(regions, '* TreatF', collapse = ' + ')))))

############################
# Testing H3a--nationalistic
############################

# No Covariates

h3a_nat_pure_att_screen<-with(basedata[basedata$Screener & (basedata$TreatG | basedata$TreatA),], lm(as.formula(paste0('Attitude ~ TreatG + ', paste(regions, '* TreatG', collapse = ' + ')))))
h3a_nat_reg_att_screen<-with(basedata[basedata$Screener & (basedata$TreatG | basedata$TreatB),], lm(as.formula(paste0('Attitude ~ TreatG + ', paste(regions, '* TreatG', collapse = ' + ')))))
h3a_nat_pure_beh_screen<-with(basedata[basedata$Screener & (basedata$TreatG | basedata$TreatA),], lm(as.formula(paste0('Behavior ~ TreatG + ', paste(regions, '* TreatG', collapse = ' + ')))))
h3a_nat_reg_beh_screen<-with(basedata[basedata$Screener & (basedata$TreatG | basedata$TreatB),], lm(as.formula(paste0('Behavior ~ TreatG + ', paste(regions, '* TreatG', collapse = ' + ')))))

# Covariates

h3a_nat_pure_att_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatG | basedata$TreatA),], lm(as.formula(paste0('Attitude ~ TreatG + ', paste(covars, '* TreatG', collapse = ' + '), '+', paste(regions, '* TreatG', collapse = ' + ')))))
h3a_nat_reg_att_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatG | basedata$TreatB),], lm(as.formula(paste0('Attitude ~ TreatG + ', paste(covars, '* TreatG', collapse = ' + '), '+', paste(regions, '* TreatG', collapse = ' + ')))))
h3a_nat_pure_beh_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatG | basedata$TreatA),], lm(as.formula(paste0('Behavior ~ TreatG + ', paste(covars, '* TreatG', collapse = ' + '), '+', paste(regions, '* TreatG', collapse = ' + ')))))
h3a_nat_reg_beh_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatG | basedata$TreatB),], lm(as.formula(paste0('Behavior ~ TreatG + ', paste(covars, '* TreatG', collapse = ' + '), '+', paste(regions, '* TreatG', collapse = ' + ')))))

############
# Testing H4
############

# No Covariates

h4_pure_att_screen<-with(basedata[basedata$Screener & (basedata$TreatH | basedata$TreatA),], lm(as.formula(paste0('Attitude ~ TreatH + ', paste(regions, '* TreatH', collapse = ' + ')))))
h4_reg_att_screen<-with(basedata[basedata$Screener & (basedata$TreatH | basedata$TreatB),], lm(as.formula(paste0('Attitude ~ TreatH + ', paste(regions, '* TreatH', collapse = ' + ')))))
h4_pure_beh_screen<-with(basedata[basedata$Screener & (basedata$TreatH | basedata$TreatA),], lm(as.formula(paste0('Behavior ~ TreatH + ', paste(regions, '* TreatH', collapse = ' + ')))))
h4_reg_beh_screen<-with(basedata[basedata$Screener & (basedata$TreatH | basedata$TreatB),], lm(as.formula(paste0('Behavior ~ TreatH + ', paste(regions, '* TreatH', collapse = ' + ')))))

# Covariates

h4_pure_att_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatH | basedata$TreatA),], lm(as.formula(paste0('Attitude ~ TreatH + ', paste(covars, '* TreatH', collapse = ' + '), '+', paste(regions, '* TreatH', collapse = ' + ')))))
h4_reg_att_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatH | basedata$TreatB),], lm(as.formula(paste0('Attitude ~ TreatH + ', paste(covars, '* TreatH', collapse = ' + '), '+', paste(regions, '* TreatH', collapse = ' + ')))))
h4_pure_beh_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatH | basedata$TreatA),], lm(as.formula(paste0('Behavior ~ TreatH + ', paste(covars, '* TreatH', collapse = ' + '), '+', paste(regions, '* TreatH', collapse = ' + ')))))
h4_reg_beh_covar_screen<-with(basedata[basedata$Screener & (basedata$TreatH | basedata$TreatB),], lm(as.formula(paste0('Behavior ~ TreatH + ', paste(covars, '* TreatH', collapse = ' + '), '+', paste(regions, '* TreatH', collapse = ' + ')))))

# Gather coefficients and confidence intervals

ci_level<-.9

ate_pure_att_screen<-sapply(list(h1_pure_att_screen, h2_pure_att_neg_screen, h2_pure_att_pos_screen, h3a_prag_pure_att_screen, h3a_ideol_pure_att_screen, h3a_nat_pure_att_screen, h4_pure_att_screen), function(x)  coef(x)[2])
ate_pure_att_covar_screen<-sapply(list(h1_pure_att_covar_screen, h2_pure_att_neg_covar_screen, h2_pure_att_pos_covar_screen, h3a_prag_pure_att_covar_screen, h3a_ideol_pure_att_covar_screen, h3a_nat_pure_att_covar_screen, h4_pure_att_covar_screen), function(x)  coef(x)[2])
ate_pure_att_covar_screen[3]<-NA # Not estimated properly given small sample
ate_reg_att_screen<-sapply(list(h1_reg_att_screen, h2_reg_att_neg_screen, h2_reg_att_pos_screen, h3a_prag_reg_att_screen, h3a_ideol_reg_att_screen, h3a_nat_reg_att_screen, h4_reg_att_screen), function(x)  coef(x)[2])
ate_reg_att_covar_screen<-sapply(list(h1_reg_att_covar_screen, h2_reg_att_neg_covar_screen, h2_reg_att_pos_covar_screen, h3a_prag_reg_att_covar_screen, h3a_ideol_reg_att_covar_screen, h3a_nat_reg_att_covar_screen, h4_reg_att_covar_screen), function(x)  coef(x)[2])
ate_reg_att_covar_screen[3]<-NA # Not estimated properly given small sample

ci_pure_att_screen<-sapply(list(h1_pure_att_screen, h2_pure_att_neg_screen, h2_pure_att_pos_screen, h3a_prag_pure_att_screen, h3a_ideol_pure_att_screen, h3a_nat_pure_att_screen, h4_pure_att_screen), function(x) confint(x, level=ci_level)[2,])
ci_pure_att_covar_screen<-sapply(list(h1_pure_att_covar_screen, h2_pure_att_neg_covar_screen, h2_pure_att_pos_covar_screen, h3a_prag_pure_att_covar_screen, h3a_ideol_pure_att_covar_screen, h3a_nat_pure_att_covar_screen, h4_pure_att_covar_screen), function(x) confint(x, level=ci_level)[2,])
ci_pure_att_covar_screen[,3]<-NA # Not estimated properly given small sample
ci_reg_att_screen<-sapply(list(h1_reg_att_screen, h2_reg_att_neg_screen, h2_reg_att_pos_screen, h3a_prag_reg_att_screen, h3a_ideol_reg_att_screen, h3a_nat_reg_att_screen, h4_reg_att_screen), function(x) confint(x, level=ci_level)[2,])
ci_reg_att_covar_screen<-sapply(list(h1_reg_att_covar_screen, h2_reg_att_neg_covar_screen, h2_reg_att_pos_covar_screen, h3a_prag_reg_att_covar_screen, h3a_ideol_reg_att_covar_screen, h3a_nat_reg_att_covar_screen, h4_reg_att_covar_screen), function(x) confint(x, level=ci_level)[2,])
ci_reg_att_covar_screen[,3]<-NA # Not estimated properly given small sample

ate_pure_beh_screen<-sapply(list(h1_pure_beh_screen, h2_pure_beh_neg_screen, h2_pure_beh_pos_screen, h3a_prag_pure_beh_screen, h3a_ideol_pure_beh_screen, h3a_nat_pure_beh_screen, h4_pure_beh_screen), function(x)  coef(x)[2])
ate_pure_beh_covar_screen<-sapply(list(h1_pure_beh_covar_screen, h2_pure_beh_neg_covar_screen, h2_pure_beh_pos_covar_screen, h3a_prag_pure_beh_covar_screen, h3a_ideol_pure_beh_covar_screen, h3a_nat_pure_beh_covar_screen, h4_pure_beh_covar_screen), function(x)  coef(x)[2])
ate_pure_beh_covar_screen[3]<-NA # Not estimated properly given small sample
ate_reg_beh_screen<-sapply(list(h1_reg_beh_screen, h2_reg_beh_neg_screen, h2_reg_beh_pos_screen, h3a_prag_reg_beh_screen, h3a_ideol_reg_beh_screen, h3a_nat_reg_beh_screen, h4_reg_beh_screen), function(x)  coef(x)[2])
ate_reg_beh_covar_screen<-sapply(list(h1_reg_beh_covar_screen, h2_reg_beh_neg_covar_screen, h2_reg_beh_pos_covar_screen, h3a_prag_reg_beh_covar_screen, h3a_ideol_reg_beh_covar_screen, h3a_nat_reg_beh_covar_screen, h4_reg_beh_covar_screen), function(x)  coef(x)[2])
ate_reg_beh_covar_screen[3]<-NA # Not estimated properly given small sample

ci_pure_beh_screen<-sapply(list(h1_pure_beh_screen, h2_pure_beh_neg_screen, h2_pure_beh_pos_screen, h3a_prag_pure_beh_screen, h3a_ideol_pure_beh_screen, h3a_nat_pure_beh_screen, h4_pure_beh_screen), function(x) confint(x, level=ci_level)[2,])
ci_pure_beh_covar_screen<-sapply(list(h1_pure_beh_covar_screen, h2_pure_beh_neg_covar_screen, h2_pure_beh_pos_covar_screen, h3a_prag_pure_beh_covar_screen, h3a_ideol_pure_beh_covar_screen, h3a_nat_pure_beh_covar_screen, h4_pure_beh_covar_screen), function(x) confint(x, level=ci_level)[2,])
ci_pure_beh_covar_screen[,3]<-NA # Not estimated properly given small sample
ci_reg_beh_screen<-sapply(list(h1_reg_beh_screen, h2_reg_beh_neg_screen, h2_reg_beh_pos_screen, h3a_prag_reg_beh_screen, h3a_ideol_reg_beh_screen, h3a_nat_reg_beh_screen, h4_reg_beh_screen), function(x) confint(x, level=ci_level)[2,])
ci_reg_beh_covar_screen<-sapply(list(h1_reg_beh_covar_screen, h2_reg_beh_neg_covar_screen, h2_reg_beh_pos_covar_screen, h3a_prag_reg_beh_covar_screen, h3a_ideol_reg_beh_covar_screen, h3a_nat_reg_beh_covar_screen, h4_reg_beh_covar_screen), function(x) confint(x, level=ci_level)[2,])
ci_reg_beh_covar_screen[,3]<-NA # Not estimated properly given small sample

# Appendix Figure 9

lwd<-2
adjust<-0.15
nvars<-length(ate_pure_att_screen)
cex<-1
cex.point<-1.5

pdf(file='combo_att_screen.pdf',width=6.5,height=8)
par(mar=c(3,3,1,1)+.01)
plot(ate_pure_att_covar_screen, nvars:1+0.5*adjust, type='p',axes=F,frame.plot=T,xlab='',ylab='', pch=19,xlim=c(-1,1), ylim=c(1-1.5*adjust,nvars+1.5*adjust),main='',sub='',cex=cex.point)
Axis(side=1,at=seq(-1, 1, 0.25),labels= seq(-1, 1, 0.25),cex.axis=cex)
Axis(side=2,at=c((nvars+1):0),labels=c('','Military\nCrimes','Base\nExpansion','Base\nReduction','Pragmatic\nFraming','Ideological\nFraming','Nationalistic\nFraming','Local\nLeaders',''),las=3,cex.axis=cex, tick=F)
abline(v=0,lwd=1,lty=2)
segments(ci_pure_att_covar_screen[1,],nvars:1+0.5*adjust, ci_pure_att_covar_screen[2,],nvars:1+0.5*adjust,lwd=lwd)
segments(ci_pure_att_screen[1,],nvars:1+1.5*adjust, ci_pure_att_screen[2,],nvars:1+1.5*adjust,lwd=lwd)
segments(ci_reg_att_screen[1,],nvars:1-0.5*adjust, ci_reg_att_screen[2,],nvars:1-0.5*adjust,lwd=lwd)
segments(ci_reg_att_covar_screen[1,],nvars:1-1.5*adjust, ci_reg_att_covar_screen[2,],nvars:1-1.5*adjust,lwd=lwd)
points(ate_pure_att_screen, nvars:1+1.5*adjust, pch=21, bg='white', cex=cex.point)
points(ate_reg_att_screen, nvars:1-0.5*adjust, pch=24, bg='white', cex=cex.point)
points(ate_reg_att_covar_screen, nvars:1-1.5*adjust, pch=17, cex=cex.point)
legend('topleft',pch=c(21,19,24,17),legend=c('Pure / No','Pure / Yes','Regular / No','Regular / Yes'), title='Control / Covariates:', inset=0.04)
dev.off()

# Appendix Figure 10

pdf(file='combo_beh_screen.pdf',width=6.5,height=8)
par(mar=c(3,3,1,1)+.01)
plot(ate_pure_beh_covar_screen, nvars:1+0.5*adjust, type='p',axes=F,frame.plot=T,xlab='',ylab='', pch=19,xlim=c(-.5,.5), ylim=c(1-1.5*adjust,nvars+1.5*adjust),main='',sub='',cex=cex.point)
Axis(side=1,at=seq(-.5, .5, 0.125),labels= seq(-.5, .5, 0.125),cex.axis=cex)
Axis(side=2,at=c((nvars+1):0),labels=c('','Military\nCrimes','Base\nExpansion','Base\nReduction','Pragmatic\nFraming','Ideological\nFraming','Nationalistic\nFraming','Local\nLeaders',''),las=3,cex.axis=cex, tick=F)
abline(v=0,lwd=1,lty=2)
segments(ci_pure_beh_covar_screen[1,],nvars:1+0.5*adjust, ci_pure_beh_covar_screen[2,],nvars:1+0.5*adjust,lwd=lwd)
segments(ci_pure_beh_screen[1,],nvars:1+1.5*adjust, ci_pure_beh_screen[2,],nvars:1+1.5*adjust,lwd=lwd)
segments(ci_reg_beh_screen[1,],nvars:1-0.5*adjust, ci_reg_beh_screen[2,],nvars:1-0.5*adjust,lwd=lwd)
segments(ci_reg_beh_covar_screen[1,],nvars:1-1.5*adjust, ci_reg_beh_covar_screen[2,],nvars:1-1.5*adjust,lwd=lwd)
points(ate_pure_beh_screen, nvars:1+1.5*adjust, pch=21, bg='white', cex=cex.point)
points(ate_reg_beh_screen, nvars:1-0.5*adjust, pch=24, bg='white', cex=cex.point)
points(ate_reg_beh_covar_screen, nvars:1-1.5*adjust, pch=17, cex=cex.point)
legend('topleft',pch=c(21,19,24,17),legend=c('Pure / No','Pure / Yes','Regular / No','Regular / Yes'), title='Control / Covariates:', inset=0.04)
dev.off()
